Bienvenid@s a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Enlaces a videos de las clases:
Documentación:
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
##
## Attaching package: 'plotly'
##
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
##
## The following object is masked from 'package:stats':
##
## filter
##
##
## The following object is masked from 'package:graphics':
##
## layout
##
##
##
## Attaching package: 'gridExtra'
##
##
## The following object is masked from 'package:dplyr':
##
## combine
Elimine eval=FALSE para ejecutar las celdas.
estimador <- function(sigma=0) {
for (n in x) {
e <- rnorm(1, 0, sigma) + (1/n) * sum(rbernoulli(cte)[1:n])
data <<- c(data, e)
}
}
mean_estimator_bias <- function(estimator_values, param) {
return(mean(estimator_values) - param)
}
# sigma = 0
cte <- rep(0.5, 1000)
x <- 1:1000
data <- c()
estimador()
## Warning: `rbernoulli()` was deprecated in purrr 1.0.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot(data, main="Estimador para Bernoulli con σ=0", xlab="n", ylab="Estimador")
abline(h=0.5, col="red")
legend('topright', "Valor Real", col="red", lty = c(1,1), bty ='n', cex=.75)
cat("σ=0\n")
## σ=0
cat("Bias:", mean_estimator_bias(data, 0.5), "\n")
## Bias: 0.0009507603
# sigma = 1
data <- c()
estimador(1)
plot(data, main="Estimador para Bernoulli con σ=1", xlab="n", ylab="Estimador")
abline(h=0.5, col="red")
legend('topright', "Valor Real", col="red", lty = c(1,1), bty ='n', cex=.75)
cat("σ=1\n")
## σ=1
cat("Bias:", mean_estimator_bias(data, 0.5), "\n")
## Bias: 0.009312616
# sigma = 2
data <- c()
estimador(2)
plot(data, main="Estimador para Bernoulli con σ=2", xlab="n", ylab="Estimador")
abline(h=0.5, col="red")
legend('topright', "Valor Real", col="red", lty = c(1,1), bty ='n', cex=.75)
cat("σ=2\n")
## σ=2
cat("Bias:", mean_estimator_bias(data, 0.5), "\n")
## Bias: -0.01850244
# sigma = 4
data <- c()
estimador(4)
plot(data, main="Estimador para Bernoulli con σ=4", xlab="n", ylab="Estimador")
abline(h=0.5, col="red")
legend('topright', "Valor Real", col="red", lty = c(1,1), bty ='n', cex=.75)
cat("σ=4\n")
## σ=4
cat("Bias:", mean_estimator_bias(data, 0.5), "\n")
## Bias: 0.07935004
Respuesta
Respuesta
Se evidencia claramente como al introducir valores aleatorios que siguen la distribución normal centrada en 0 pero con un \(\sigma > 0\) el estimador pierde consistencia, dado que pese a cumplir la primera condición, es decir, tener un sesgo cercano a 0, a medida que aumenta el valor de \(\sigma\) las estimaciones no tienden al valor real, lo cual es la segunda condición para considerarse consistente.
En este caso se demostró empíricamente que los estimadores con menos sesgo no necesariamente son los más precisos dado el argumento previo.
En las preguntas 2, 3 y 4 deberá trabajar con el dataset
diabetes_prediction_dataset.
datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)
Una forma de estimar la distribución de la media de una población es
a través de la realización de la sampling distribution de
una distribución cualquiera. El valor obtenido en cada sampling
distribution nos entrega un estimador de la media que posee una
determinada distribución de la población. Sabiendo esto, calcule la
sampling distribution de las columnas bmi y
blood_glucose_level, obteniendo el intervalo de confianza
de \(95\%\) para cada una de las medias
obtenidas desde la distribución. Para realizar este experimento
considere los siguientes puntos:
sampling distribution con un tamaño de muestra igual a
\(100\). Repita la obtención de la
media un número elevado de veces (recomendación \(5000\) veces).bmi y blood_glucose_level y señale que
diferencias y similitudes encuentra entre estas (en cuanto a la calidad
de la estimación, no los valores obtenidos). De una explicación de por
qué cree que se dan las similitudes/diferencias.Hints:
bmi y blood_glucose_level a una misma
escala.sampling distribution podría serle
útil el comando sample().Respuesta:
# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 10000
# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_izq_bmi = vector()
list_interval_der_bmi = vector()
list_typeCI_bmi = vector()
prop_bmi = 0
sucess_bmi = 0
list_mean_bgl = vector()
list_interval_izq_bgl = vector()
list_interval_der_bgl = vector()
list_typeCI_bgl = vector()
prop_bgl = 0
sucess_bgl = 0
#normalizamos columnas
bmi_col <- na.omit((datos$bmi - min(datos$bmi))/(max(datos$bmi)-min(datos$bmi)))
bgl_col <- na.omit((datos$blood_glucose_level - min(datos$blood_glucose_level))/(max(datos$blood_glucose_level)-min(datos$blood_glucose_level)))
# Obtenemos la media y desviación estándar de cada columna
real_mean_bmi <- mean(bmi_col, na.rm = TRUE)
real_mean_bgl <- mean(bgl_col, na.rm = TRUE)
real_std_bmi <- sd(bmi_col, na.rm = TRUE)
real_std_bgl <- sd(bgl_col, na.rm = TRUE)
# Sampling distribution, calculo del intervalo de confianza y proporción.
for (i in 1:n_muestras) {
#hacemos los samples
sample_bmi <-sample(bmi_col, tamano_muestra,replace=TRUE)
sample_bgl <- sample(bgl_col, tamano_muestra,replace=TRUE)
#calculamos y guardamos las medias de cada sample
mean_sample_bmi <- mean(sample_bmi)
list_mean_bmi <- c(list_mean_bmi, mean_sample_bmi)
mean_sample_bgl <- mean(sample_bgl)
list_mean_bgl <- c(list_mean_bgl, mean_sample_bgl)
#para la fórmula de intervalo de confianza
zmedio <- qt(1 - 0.05 / 2, df =tamano_muestra-1)
lim_izq_bmi <- mean_sample_bmi - (zmedio*real_std_bmi/sqrt(tamano_muestra))
lim_der_bmi <- mean_sample_bmi + (zmedio*real_std_bmi/sqrt(tamano_muestra))
list_interval_izq_bmi <- c(list_interval_izq_bmi, lim_izq_bmi)
list_interval_der_bmi <- c(list_interval_der_bmi, lim_der_bmi)
if(lim_izq_bmi <= real_mean_bmi & lim_der_bmi >= real_mean_bmi){
sucess_bmi <- sucess_bmi+1
list_typeCI_bmi <- c(list_typeCI_bmi,TRUE)
}else{
list_typeCI_bmi <- c(list_typeCI_bmi,FALSE)
}
lim_izq_bgl <- mean_sample_bgl - (zmedio*real_std_bgl/sqrt(tamano_muestra))
lim_der_bgl <- mean_sample_bgl + (zmedio*real_std_bgl/sqrt(tamano_muestra))
list_interval_izq_bgl <- c(list_interval_izq_bgl, lim_izq_bgl)
list_interval_der_bgl <- c(list_interval_der_bgl, lim_der_bgl)
if(lim_izq_bgl <= real_mean_bgl & lim_der_bgl >= real_mean_bgl){
sucess_bgl <- sucess_bgl+1
list_typeCI_bgl <- c(list_typeCI_bgl,TRUE)
}else{
list_typeCI_bgl <- c(list_typeCI_bgl,FALSE)
}
}
# Generamos dataframes para plotear mas facilmente los datos.
dfp2_bmi <- data.frame(y = seq(1, 100),Mean_bmi=list_mean_bmi[1:100], Izq_CI_bmi = list_interval_izq_bmi[1:100], Der_CI_bmi = list_interval_der_bmi[1:100], In_Out_bmi=list_typeCI_bmi[1:100])
dfp2_bgl <- data.frame(y = seq(1,100), Mean_bgl=list_mean_bgl[1:100], Izq_CI_bgl = list_interval_izq_bgl[1:100], Der_CI_bgl = list_interval_der_bgl[1:100], In_Out_bgl=list_typeCI_bgl[1:100])
#calculamos las proporciones
prop_bmi<- sucess_bmi/n_muestras
prop_bgl <- sucess_bgl/n_muestras
# Plot de Intervalos de confianza
mu_bmi <- mean(dfp2_bmi$Mean_bmi)
ggplot(dfp2_bmi, aes(x = Mean_bmi, y =y)) +
geom_errorbarh(aes(xmin = Izq_CI_bmi, xmax = Der_CI_bmi , color = factor(In_Out_bmi)), height = 0.2) +
geom_point() +
geom_vline(xintercept = mu_bmi, color = "black", linetype = "solid") + # Línea en la media de BMI
labs(title = "Forest Plot para BMI", x = "BMI", y = "Muestras") +
theme_minimal()
# Plot de Intervalos de confianza
mu_bgl <- mean(dfp2_bgl$Mean_bgl)
ggplot(dfp2_bgl, aes(x = Mean_bgl, y = y)) +
geom_errorbarh(aes(xmin = Izq_CI_bgl, xmax = Der_CI_bgl , color = factor(In_Out_bgl)), height = 0.2) +
geom_point() +
geom_vline(xintercept = mu_bgl, color = "black", linetype = "solid") + # Línea en la media de BMI
labs(title = "Forest Plot para BGL", x = "BGL", y = "Muestras") +
theme_minimal()
barplot(c(sucess_bmi, n_muestras - sucess_bmi),
names.arg = c("Exitos", "Fracasos"),
col = c("skyblue", "pink"),
main = "BMI",
ylab = "cantidad")
barplot(c(sucess_bgl, n_muestras - sucess_bgl),
names.arg = c("Exitos", "Fracasos"),
col = c("skyblue", "pink"),
main = "BGL",
ylab = "cantidad")
# Plot de proporción de CI
cant.unos <- cumsum(list_typeCI_bgl)
cant.unos2 <- cumsum(list_typeCI_bmi)
# Después calculamos la proporción de 1s en cada punto
# (normalizamos)
prop.unos <- cant.unos/seq(from=1, to=length(cant.unos), by=1)
prop.unos2 <- cant.unos2/seq(from=1, to=length(cant.unos2), by=1)
# Gráfico de la convergencia
plot(prop.unos, xlab = "Cantidad de samples", ylab = "Probabilidad", type = "l", col = "pink", lwd = 2, ylim = c(0, 1))
lines(prop.unos2, col = "purple", lwd = 2)
abline(h = 0.95, col = "red", lty = 2)
legend("bottomright", legend = c("BGL", "BMI"), col = c("pink", "purple"), lwd = 2)
Respuesta
En cuanto a la proporción de los datos se evidencia que estos concuerdan con la confianza del 95%, mientras que respecto a la calidad de la estimación, se puede apreciar visualmente que para el caso de bgl los valores varían un poco más. Para compararlos calcularemos el error estándar: \(SE=\frac{\sigma}{\sqrt{n}}\). En este caso n sería para ambos igual a 100. Para bmi tenemos: \(\frac{0.07746}{10}\) Para bgl tenemos: \(\frac{0.1850}{10}\)
Es decir, SEbmi < SEbgl, lo cual explica la observación previamente mencionada. Las estimaciones de la media para bmi son más precisas que las de bgl, ya que el error estándar de bmi es menor. Un error estándar más bajo indica que las medias muestrales están más cercanas a la media real de la población. Esto implica intervalos de confianza más estrechos y más precisión en la estimación.
Siguiendo lo anterior mencionado, efectivamente se puede observar que los largos de los intervalos en el caso de bgl son un poco mayores que los de bmi, esto también se explica por el hecho de que stdbgl>stdbmi. Al tener los datos una mayor variabilidad afecta a los resultados haciendo que se dificulte un poco más estimar de forma más acertada.
El objetivo de esta parte será estimar la media bmi
mediante MLE. Para responder la pregunta realice los siguientes
puntos:
bmi y obtenga la media y varianza de los datos.\[ loglike(\mu) = -n\log\sigma -\frac{nS^2}{2\sigma^2} -\frac{n(\bar{X}-\mu)^2}{2\sigma^2} \] \[ S = n^{-1}\sum^n_{i=1} X_i - \bar{X}^2 \]
persp() o
filled.contour()). Considere el histograma de
bmi y el valor de \(\hat{\mu}\) para definir un rango adecuado
para la media, y que \(\sigma\) está
entre 1 y 7.nlminb() sobre la función likelihood (o log-likelihood)
para encontrar el óptimo de \(\mu\) y
\(\sigma\) computacionalemte.nlminb() con el resultado teórico y la media muestral, y
comente sus resultados.Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).
Respuesta
# Histograma
#hist(datos$bmi,main = "Histograma de BMI",breaks=20, xlab = "BMI", ylab = "Frecuencia", col = "pink", border = "black")
ggplot(datos, aes(x = bmi)) +
geom_histogram(binwidth = 1, fill = "pink", color = "black", alpha = 0.7) +
labs(title = "Histograma de BMI", x = "BMI", y = "Frecuencia") +
theme_minimal()
media_p3<-mean(datos$bmi)
varianza_p3 <- var(datos$bmi)
Cálculo de valor teórico del estimador de máxima verosimilitud \(\hat{\mu}\) para la media, asumiendo una distribución normal y fijando \(\sigma^2\) en la varianza muestral. Derivamos con respecto a mu. \[ \frac {\partial Loglike}{\partial \mu}= \frac{\partial}{\partial \mu}(-\frac{n*(\bar{X}^{2}-2\bar{X}\mu+\mu^2)}{2\sigma^2}) = \frac{n*(\bar{X} - \mu)}{\sigma^2} \] Y ahora procedemos a igualar a cero:
\[ \frac{n*(\bar{x} - \mu)}{\sigma^2}=0 \Rightarrow \bar{X} - \mu=0 \Rightarrow \bar{X} = \mu \]
Así se concluye que el valor teórico del estimador de máxima verosimilitud \(\hat{\mu}\) para la media es \(\bar{X} = \mu\)
# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
n<- length(datos$bmi)
X_raya <- mean(datos$bmi)
S_cuadrado <- sum((datos$bmi - X_raya)^2) / n
resultado <- (-n * log(sigma)) - ((n * S_cuadrado )/ (2 * sigma^2)) -
(n * ((X_raya - mu)^2) / (2 * (sigma^2)))
return (-resultado)
}
ll_plot <- Vectorize(ll_plot)
# definir espacio donde se va a evaluar ll_plot
mu_m <- seq(27,28,by=0.1)
sigma_m <- seq(1,7, by=0.1)
likelihood_m <- outer(X=mu_m, Y=sigma_m, ll_plot)
filled.contour(x=mu_m, y=sigma_m, z=likelihood_m, xlab=expression(mu),
ylab=expression(sigma), color = topo.colors)
# definir espacio donde se va a evaluar ll_plot
mu_m <- seq(27,28,by=0.1)
sigma_m <- seq(6,7, by=0.1)
likelihood_m <- outer(X=mu_m, Y=sigma_m, ll_plot)
filled.contour(x=mu_m, y=sigma_m, z=likelihood_m, xlab=expression(mu),
ylab=expression(sigma), color = topo.colors)
likelihood_mean <- function(param) {
# Definimos los parametros de entrada de la funcion
mu <- param[1]
sigma <- param[2]
# Definimos la likelihood como la suma logaritmica de la función de densidad
ll_plot(mu,sigma)
# ...
}
# Agregue el rango donde operará nlminb
nlminb(objective=likelihood_mean, start=c(27,6), lower=c(27, 1), upper=c(28, 7))
## $par
## [1] 27.320776 6.636744
##
## $objective
## [1] 239262.2
##
## $convergence
## [1] 0
##
## $iterations
## [1] 10
##
## $evaluations
## function gradient
## 19 28
##
## $message
## [1] "relative convergence (4)"
Al comparar los valores de la media real y la teórica se tiene que son las misma, además, la obtenida con max likelihood nos dió 27.32077. Por lo que se establece que la estimación es muy buena.
Con respecto a los resultados obtenidos, vemos que el gráfico generado centra el mínimo cerca del mínimo real, lo cual indica que efectivamente está bien hecho y es el resultado que se esperaba.
En esta pregunta se busca el error estándar y el intervalo de
confianza para la varianza de la columna HbA1c_level.
Las actividades por realizar son las siguientes:
Respuesta:
# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()
for (it in 1:B) {
output <- c(output, var(sample(datos$HbA1c_level, largo, replace=T)))
}
# Gráfico de puntos
se <- sd(output)
valor_real <- mean(output)
plot(output, main="Varianzas obtenidas con Bootstrap", xlab="iteración", ylab="varianza")
abline(h=valor_real, col="red")
legend('topright', "Valor Real", col="red", lty = c(1,1), bty ='n', cex=.75)
# Histograma
h <- hist(output, breaks=20, plot=F)
h$density <- h$counts / sum(h$counts)
plot(h, freq=F, main="Histograma del muestreo con Bootstrap", xlab="varianza", ylab="frecuencia", col="pink")
x_vals <- seq(min(output), max(output), length = 1000)
y_vals <- dnorm(x_vals, mean(output), sd(output))
y_vals_scaled <- y_vals * sum(h$density) * diff(h$mids[1:2])
lines(x_vals, y_vals_scaled, col = "blue", lwd=2)
abline(v=valor_real, col="red")
legend('topright', c("Valor Real", "Sample dnorm"), col=c("red", "blue"), lty = c(1, 1), cex=.75)
# Obtenemos error e intervalos de confianza
se_vars <- se
alpha <- 0.05
# límite inferior
l.CI <- quantile(output, alpha/2)
# límite superior
u.CI <- quantile(output, 1-alpha/2)
sprintf('El intervalo de confianza del 95%% de las varianzas es: (%.5f, %.5f)', l.CI , u.CI)
## [1] "El intervalo de confianza del 95% de las varianzas es: (1.13580, 1.15700)"
sprintf('El SE de las varianzas es: (%.5f)', se_vars)
## [1] "El SE de las varianzas es: (0.00542)"
Respuesta
Realizar Bootstrap sobre la columna de datos HbA1c_level
genera un histograma tal que se aproxima a una distribución normal,
mostrando además que el error estándar de la varianza es muy bajo, con
lo cual se puede inferir que el estimador de la varianza de datos tiene
una buena precisión y es posible que los datos originales tengan poca
variabilidad.
El siguiente código genera una regresión lineal de bmi
con respecto a age, es decir, una función lineal de
age, \(y(age) = b + m\cdot
age\), que pretende determinar el valor de bmi.
linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)
##
## Call:
## lm(formula = bmi ~ age, data = datos)
##
## Coefficients:
## (Intercept) age
## 23.15536 0.09945
m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]
ggplot() +
geom_point(aes(x=datos$age, y=datos$bmi)) +
geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
ggtitle("Regresión lineal") +
ylab("bmi") +
xlab("Age")
Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap? ¿Un valor bajo en el error significa que la regresión es buena?
# Bootstrap
B <- 500
largo <- 1000
output_m <- c()
output_b <- c()
for (it in 1:B) {
output_m <- c(output_m, lm(bmi ~ age, datos[sample(nrow(datos), largo, replace=T), ])$coefficients["age"])
output_b <- c(output_b, lm(bmi ~ age, datos[sample(nrow(datos), largo, replace=T), ])$coefficients["(Intercept)"])
}
# Obtenemos error e intervalos de confianza
se_m <- sd(output_m)
alpha <- 0.05
# límite inferior
l.CI <- quantile(output_m, alpha/2)
# límite superior
u.CI <- quantile(output_m, 1-alpha/2)
sprintf('El intervalo de confianza del 95%% de m es: (%.3f, %.3f)', l.CI , u.CI)
## [1] "El intervalo de confianza del 95% de m es: (0.082, 0.116)"
sprintf('El SE de m es: (%.5f)', se_m)
## [1] "El SE de m es: (0.00841)"
# Obtenemos error e intervalos de confianza
se_b <- sd(output_b)
alpha <- 0.05
# límite inferior
l.CI <- quantile(output_b, alpha/2)
# límite superior
u.CI <- quantile(output_b, 1-alpha/2)
sprintf('El intervalo de confianza del 95%% de b es: (%.3f, %.3f)', l.CI , u.CI)
## [1] "El intervalo de confianza del 95% de b es: (22.365, 23.922)"
sprintf('El SE de b es: (%.3f)', se_b)
## [1] "El SE de b es: (0.410)"
Respuesta
Es posible evidenciar que el error estándar para ambas variables es
relativamente pequeño, indicando que al calcular los parámetros
m y b según las muestras de Bootstrap, estos
resultan en parámetros similares para la regresión lineal según el
estimador, esto, dado que al tomar samples de tamaño mil para cada
iteración, los datos en general tomarán una distribución similar a los
originales pese al muestreo realizado, resultando efectivamente en una
buena estimación de estos valores en cada iteración. Que este error sea
bajo solamente indica que las estimaciones de m y
b son consistentes y estables, pero no necesariamente que
la regresión lineal sea “buena” en el sentido de capturar adecuadamente
la relación entre estas variables, para verificar de mejor forma si esta
regresión es buena sería útil calcular el coeficiente de determinación
\(R^2\) y analizar su valor.
A work by CC6104